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I. Introduction 

Mathematical modeling of the interactions of large, high- 
voltage solar arrays with the space plasma in Low Earth Orbit (LEO) 
has been conducted by the P.I. since 1984. An early series of 
models sought to identify the main features of the electrical 
response of a model solar array/power system to arcing at points of 
negative bias on the array relative to the plasma potential (refs. 1- 

4) . In addition to characterizing array responses to various arcs, 
this early work suggested that a description of the electrically 
floating state of an array prior to an arc requires a correct 
description of the plasma sheath that surrounds the array. 

One approach to the treatment of the plasma sheath near a 
solar array is contained in the two-fluid, warm plasma model (ref. 

5) . If magnetic fields are ignored, the plasma equations in this 
model are the equations of continuity and of motion for the 
electrons and the ions and Poisson's Equation connecting the 
electron and ion charge densities with the electrostatic potential. In 
three dimensions, for irrotational flows, these comprise five partial 
differential equations in five scalar fields: the electron and ion 
charge densities, the electron and ion velocity potentials and the 
electrostatic potential. 

An active solar array has a distribution of voltage bias 
relative to plasma ground over its surface. Solar cell layouts 
designed to minimize vxB forces necessitate that some points of 
sharp voltage discontinuity will be present on the array. Points of 
shape discontinuity, such as edges and comers, also occur. Proper 
description of the sheath near such discontinuities can be expected 


to require that the sheath equations be solved in their multi- 
dimensional form. Approximations based on one-dimensional 
assumptions embedded in a sheath model are potential sources of 
error. 


In the limiting cases when surface bias voltages are either 
low or high as compared with the electron and ion thermal energies 
divided by the electron charge, the equations of the two-fluid, 
warm plasma model yield the Debye (low bias) or Child-Langmuir 
(high bias) approximations, respectively. A fully three-dimensional 
treatment of the sheath equations for arbitrarily shaped and biased 
objects is desirable but is difficult to achieve. However, some 
insight to sheath structure can be gained by less ambitious 
calculations. Herein we present calculations in two dimensions, for 
special geometries and in the limiting, Debye and Child-Langmuir 
approximations. 

II. Previous Work 

In the summer of 1987 the P.I. began a study of the electron 
sheath near surfaces of a model high-voltage solar array (ref. 6). 

An evaluation was conducted of the NASCAP/LEO computer model. 
This is a three-dimensional model designed to calculate electrostatic 
potential distributions, particle fluxes and floating conditions near 
objects of chosen shape, composition and voltage biasing in LEO 
(refs. 7-8). The model gives an extensive picture of the plasma 
sheaths near voltage-biased structures. 

Calculational features, possibly designed to make a three 
dimensional model tractable, limit NASCAP LEO’s descriptive power 
somewhat. For some kinds of objects the model has rather low 
resolution in crucial places. The biased object under study is 
confined to a 17x17x33 pt. inner grid. The next outer grid has half 
the resolution of the inner grid. Thus, to obtain the highest 
resolution available for the sheath immediately above object 
surfaces, the object must be kept well within the inner grid volume. 
Few grid points can be deployed near fine surface features of the 
object. 


NAS CAP/LEO contains an essentially one-dimensional 
treatment of the relation between the electrostatic potential, 
velocity and charge density. This permits the electrostatic 
potential to be calculated from Poisson’s Equation directly, which is 
a great convenience. Particle fluxes are calculated using particle- 
tracking once the electrostatic potential has been calculated. This 
treatment is likely to be deficient near multi-dimensional features 
such as discontinuities and does not take particle-particle 
interactions into account in the particle flux calculations. 

A plasma fluid model should be free of the above limitations 
provided that it's equations are solved in their true multi- 
dimensional form and symmetries are imposed that permit use of a 
higher resolution grid. To test this, four plasma sheath problems 
were defined for study by the P.I. in 1987. In each of two 
geometries (Fig. 1), an equipotential, rectangular bar and a sheet of 
alternating-bias strips, and in each of the two limiting 
approximations, Debye and Child-Langmuir, the equations of the 
warm plasma model were developed and programmed for 
computer solution in two dimensions. During the summer of 1988 
preliminary solutions of the Debye equations were obtained (ref. 9). 
During the summer of 1989 final Debye solutions were obtained, 
comparisons to NASCAP/LEO results made and work on the Child- 
Langmuir cases begun. 
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GEOMETRIES: Fig. l.a is a cross-section taken perpendicular to a 

plane containing strips of voltages alternately Vq volts above and 



below V 0 o- The strips run from plus to minus infinity in the z 
direction and, side by side, to plus and minus infinity in the x 
direction. Shown is a region half a strip wide on either side of the 
origin. Fig. l.b shows one quadrant of a cross-section through a 
rectangular bar at voltage V 0 o- The bar extends to plus and minus 
infinity in the z direction. 


III. Work Conducted in 1989 under NAG 3-576 
A. Debye Approximation 
1) Sheath Equations 

Of the two limiting approximations of the fluid model, the 
Debye Approximation contains the simpler field equations. The 
Equation of Motion can be integrated exactly and yields scalar, 
exponential relations between the particle densities and the 
electrostatic potential. These relations are then used in Poisson’s 
Equation to produce a single equation for the electrostatic potential. 
Once Poisson's Equation is solved, the particle densities can be used 
in the Equation of Continuity to determine the velocity field of the 
fluid flow. 

The field equations of the Debye Approximation (eO«kT) as 
derived in the warm plasma model (ref. 9) are as follows: 

V 2 <X> - <D/(Xci) 2 = 0 (Poisson's Equation) 

V-(NV) = 0 (Equation of Continuity) (1) 

VN/N = -qVO /kT (Equation of Motion) 

where O is the electrostatic potential, A,<j is the Debye Length, N is 
the number density of a species (electrons or ions), q is the charge 
of the species particles and V is the velocity field associated with a 
species. The solution of the third of these equations is Ni on = 
N 0 exp(-e<I>/kT) and N e iectron = N 0 exp(ed>/kT), where e is the 
magnitude of the charge on the electron. Using e<D/kT«l, for the 
electrons and Nj 0n = N 0 for the ions, the above form of Poisson's 
Equation results. For the ions, the assumed zero velocity field 


trivially solves the Equation of Continuity. For the electrons, with 
the assumption of irrotational flow, i.e., VxV = 0, the Equation of 
Continuity becomes: V^F+V(eO/kT)-VF = 0 where VF = V. 

These results may be put in dimensionless form with the 
substitutions: X.dV = V', e<I>/kT = <J>' and n = N/N 0 , where N 0 is the 
species density in the ambient plasma. The velocity units are 
arbitrary and the velocity V' is derived from VF = V'. With these 
definitions the field equations become: 

(V’)2<D’ = <X>' (Poisson's Equation) 

V^f’+V'O’ V'F = 0 (Equation of Continuity for Electrons) (2) 
n = exp(O') (Relative Electron Density) 

2) Numerical Solutions 

A two-dimensional, numerical model achieving solutions of 
these equations using Central-Difference methods was developed in 
the summers of 1987 and 1988. Electrostatic potential 
distributions and electron flow fields have been modeled for both 
the edge and the strip geometries. In the case of the strip 
geometry, an analytical solution to Poisson's Equation was obtained. 
During the summer of 1989 these solutions were fully implemented 
and compared with output from NASCAP/LEO. 

To accomplish this comparison, NASCAP/LEO was run using 
three-dimensional objects spanning the full length of the 33 point 
axis of the inner grid. The objects were a long rectangular bar and 
a set of adjacent, plane rectangles. Cross-sections through the 
centers of these objects provided two-dimensional data against 
which runs of our model were compared. 

A computer program written by a student assistant, Bishwa 
Basnet, was used to read the NASCAP/LEO output files and put 
them in a form suitable for contour plotting with S, which is the 
graphics package we use to plot the output from the fluid model. 
Thus, contour plots of our output and the NASCAP/LEO output were 


made directly comparable. Results of the NASCAP/LEO electrostatic 
potential calculations were generally in agreement with fluid model 
results except in regions of the sheaths near the corner of the edge 
in the edge geometry or near the surface voltage discontinuity in 
the strip geometry. In such regions, our calculations were able to 
yield somewhat more detail than those of NASCAP/LEO. 

3) NASCAP/LEO Code Comparisons 

Figures 2-5 show comparisons in the Debye Approximation 
between NASCAP/LEO electrostatic potentials and electrostatic 
potentials calculated using the fluid model. In these figures 
voltages have been normalized to kT. Plasma parameters used to 
calculate these results are as follows: 

Electron Temperature = 0.1 e.v. 

Electron Density = 10 6 /cc 

Debye Length = 0.00235 meters 

Fig. 2 is a plot of NASCAP/LEO output based on a cross section 
through the center of a rectangular bar of dimensions 4 Xd x 2Xd x 
IIXd. The width of the bar spans 12 of the available 16 grid units 
in the inner grid of the NASCAP/LEO, the height spans 6 grid units 
and the length spans all 32 of the available units. Thus the upper 
surface of the comer shown in Fig. 2 is spanned by 7 grid points 
and the right side is spanned by 4 grid points. The solution has 
been set to zero on the boundary of the first outer grid. 

Fig. 3 shows a solution of the NASCAP/LEO-type Poisson's 
Equation (eq. 9) using the fluid model. Our calculations reflect a 
grid of 21 x 21 points and 3 Xd = one grid space. Thus our grid is 
identical to that used by NASCAP/LEO in Fig. 2. One can see that, 
except for the placement of the zero boundary (our zero is on the 
upper and right hand extreme boundaries), the equipotential lines 
in the two figures fall almost exactly on top of each other. There is 
excellent agreement between our calculations and NASCAP/LEO's in 
these figures. 

Fig. 4 is a plot of NASCAP/LEO output based on four strips 
occupying the median plane of the inner grid. Each strip is 8 grid 


spaces wide and 16 grid spaces long. The voltages of the strips 
alternate, beginning with a low voltage strip and ending with a high 
voltage strip. The arrangement is shown below in Fig. 6. Fig. 4 
shows a vertical cross-section centered on the discontinuity 
between the first low-voltage strip and the first high-voltage strip 
and perpendicular to the line of discontinuity. Only 9 grid points 
span the x axis of the calculation space, which is shown by the box. 
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Fig. 5 shows a calculation based on the same grid adopted in 
Fig. 4 but using the analytical solution of Poisson’s Equation for the 
electrostatic potential in the strip geometry (ref. 6). The two 
contour plots are similar. They differ in the placement of the 
potential zero, which is at a numerically determined distance above 
the surface in Fig. 4 and is at infinity for the calculation shown in 
Fig. 5. This artificially compresses the equipotentials toward the 
surface in Fig. 4 relative to those in Fig.5. 

The lack of good resolution in both figures leads to some 
differences near the voltage discontinuity. Poor resolution is forced 
by the the need to have four adjacent strips in the inner grid of 
NASCAP/LEO's calculation space. We were unable to achieve 
convergent runs of NASCAP/LEO using only two strips. Even using 
the longest available dimension of the inner grid, only 8 grid units 
per strip could be accommodated. The analytical solution shown in 
Fig.5 is, in principle, exact but was calculated and plotted on the 
sparse grid of NASCAP/LEO for comparison purposes. Thus, Fig. 5 
has the same lack of resolution as does Fig. 4. 

The results plotted in Figures 2-5 suggest that, in the Debye 
Approximation, our numerical approach replicates that of 


NASCAP/LEO within the precision of that model and for the 
geometries we have chosen to analyze. Where resolution is 
sufficient in both approaches, agreement is good. Where 
differences occur, poor resolution prohibits taking the differences 
too seriously. 

4) Fluid Model Runs 

Calculations in both the strip geometry and in the edge 
geometry using the fluid model were fully reported in our March, 
1989 Interim Report (ref. 9). 


B. Child-Langmuir Approximation 
1) Sheath Equations 

The Child-Langmuir Approximation describes sheath regions 
where e<I> » kT. Since kT « 1 e.v. in LEO this approximation is the 
appropriate one for modeling high-voltage solar arrays deployed 
there. Although the Equation of Motion (eqs. 2) can be integrated 
exactly, as in the Debye Approximation, the result is a quadratic 
relationship between the magnitude of the fluid velocity and the 
electrostatic potential. This is the relationship used in 
NASCAP/LEO. However, no scalar relationship between the 
electrostatic potential and the charge density arises in more than 
one dimension, as is assumed in NASCAP/LEO. Thus, in the fluid 
model, Poisson's Equation does not separate from the other 
equations as an equation for the electrostatic potential alone. The 
equations of motion and continuity remain coupled to Poisson's 
Equation and must be solved simultaneously with Poisson's 
Equation. 

If attention is focused mainly on that part of the sheath 
nearest the high voltage surfaces, where ions are excluded, the ion 
contribution to the sheath equations can be ignored and the sheath 
becomes an electron sheath. The field equations for the electron 
sheath in the Child-Langmuir Approximation, as derived in the 
warm plasma model (ref. 6) and assuming irrotational flows, are as 
follows: 
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V 2 <I> - eN/e 0 = 0 (Poisson’s Equation) 

V-(NV) = 0 (Equation of Continuity) (3) 

mV 2 /2 - e<D =0 (Equation of Motion - Energy Eq.) 

where N is the electron density, V is the electron velocity field and 
O is the electrostatic potential. The ion density differs from zero 
only on the outer boundary region of the sheath where the charge 
density term in Poisson's Equation must reflect the presence of the 
ions. Notice that the Equation of Motion also expresses conservation 
of energy for single particles. 

The appropriate scale length in the high-voltage sheath is the 
Child-Langmuir length, Dcl. given by: 

D C L 2 = (4e 0 /9)(47r/ekT)i/2<I)o 3/2 /No (4.a) 

In comparison to the Debye length, Xd: 

D C l 2 Ad 2 = (4/9)(47t)i/2(eO 0 /kT)3/2 (4.b) 

where O 0 is the surface potential, i.e., the value of the voltage on 
the surface of the edge or the average voltage of the alternating 
strips. One can see that, if ed> 0 »kT, the Child-Langmuir length 
greatly exceeds the Debye length. Thus, because the typical high- 
voltage sheath must be calculated out from the surfaces at least one 
Child-Langmuir length and only a finite number of calculational 
grid points can be deployed in such a calculation, features of the 
sheath or the surface having the scale of the Debye length are 
unlikely to be resolved. A sheath thickness of the order of Dcl is 
suggested by the result that Dcl is the thickness of the one- 
dimensional sheath above a uniformly biased, plane surface. 

As in the Debye Approximation, dimensionless forms of the 
sheath equations can be found. However, the appropriate norm 
voltage is not kT/e but O 0 and the appropriate length scale is not Xd 
but Dcl- The transformation to dimensionless units is accomplished 


by making the substitutions: Dcl v = V', O’ = O/O 0 , V’ = V/V 0 and n 

= N/N 0 , where V Q = (2eO 0 /m), m is the mass of an electron and V' is 
derived from V' = V'F. With these definitions, and using eq.(4.b), 
the field equations become: 

(V’)2<d’ - (4/9y)(kT/eO 0 )n = 0 (Poisson’s Equation) 

nV' 2 F’+V'n-V'F = 0 (Equation of Continuity) (5) 

(V*) 2 = (Equation of Motion - Energy Eq.) 

where y= (kT/eO 0 ) 3/2 /(4jt) 1/2 = 4 Xd 2 /9Dcl 2 (6) 

{Note: A numerically more useful form of the Equation of Continuity 
is V 2( g F) + (g-2)V'g-V’F - FV'2g = 0, where g s ln(n).} 

2) Numerical Solutions 

In the summer of 1989, computer programs were written to 
solve the Child-Langmuir equations in two dimensions in the edge 
and strip geometries. A self-consistent field calculation for the 
electron density, electrostatic potential and electron velocity field of 
the electron sheath was planned. 

In the self-consistent approach a trial electrostatic potential is 
chosen. The relative electron density is then deduced from 
Poisson's Equation and used in the Equation of Continuity to 
produce the electron velocity potential. The magnitude of the 
velocity is calculated throughout the grid from the velocity 
potential and the Energy Equation is used to produce a new 
electrostatic potential. The new electrostatic potential and the trial 
electrostatic potential are compared. If sufficiently different, they 
are mixed and introduced as a new trial potential beginning a new 
calculation cycle. 

To produce a trial electrostatic potential, we have chosen to 
solve the non-linear Poisson's Equation that arises using 
NASCAP/LEO’s relation connecting the charge density, velocity field 


and the electrostatic potential in the sheath (ref. 8) . As given by 
Katz, Mandell and Cooke in unnormalized units, this relationship is: 

p = -(e 0 OA D 2 )[1 + (4ic) 1 / 2 (eO/kT)3/2]-i (7) 

When eO/kT « kT, this expression reduces to the correct form of 
the net charge density in the Debye Approximation, namely, p = - 
N 0 e 2 O/kT. When eO/kT » kT, the expression gives the electron 
charge density in the one-dimensional, Child-Langmuir 
Approximation, namely, p = -N 0 e/(47ieO/kT) 1 / 2 . 

Whereas the Debye limit of eq.(7) can easily be shown to be 
correct in more than one dimension, this is not true of the Child- 
Langmuir limit. As pointed out above, there is no reason to expect 
that a scalar relationship between p and exists in more than one 
dimension for this limit. Thus, while we shall use the NAS CAP/LEO 
relationship to develop the trial potentials in our Child-Langmuir 
analysis, we expect that the self-consistent procedure will produce 
different final potentials. Differences should be most evident near 
intrinsically multi-dimensional regions of the sheath such as 
surface shape and voltage discontinuities. 

If eq.(7) is written in the dimensionless units adopted in 
eqs.(5), the result is an expression for the relative electron density: 

n = Y<D’(eOo/kT)/[Y + (0’) 3/2 ] (8) 

Using this result, Poisson's Equation becomes: 

(V’)2<d' = (4<D'/9)/[y + 0'3/ 2 ] (9) 

To solve this equation numerically, we linearize the right- 
hand side, viz: 

f(u) = u/[y+u 3 / 2 ] » 3u 0 5/2 /2(y+u 0 3/ 2 ) 2 + {(y-u<> 3/2 /2)/(y +u o 3/2 ) 2 }ii (10) 
Using this expansion, we solve (V’) 2 <1>' = (4/9)f(d>‘) iteratively. 

Initially the function O' 0 is the solution of Laplace’s Equation 



(f = 0) using the chosen boundary conditions. After one iteration 
the solution, <&', is compared with the starting function, 0’ 0 . If 
sufficiently different they are mixed and the mixture introduced as 
O’o in a new iteration. This process continues until agreement 
between O’ and O' 0 is obtained to within a chosen precision. 

As Katz, et al, have found (ref. 7), to achieve numerical 
stability it is necessary to replace the value of y in eq.(9) with a 
larger value. Since y = 4 Xd 2 /9Dcl 2 » this has the effect of replacing Xd 
as the intrinsic scale length of the sheath in regions where eO<<kT. 
Thus, in such regions the model becomes unable to resolve features 
on the scale of Xd- As Katz, et al, point out, however, the grid 
spacing for a sheath surrounding a high voltage object is typically 
chosen to be many times larger than a Debye length and there is 
little point in trying to model features smaller than a grid spacing . 

A second potential instability arises when the coefficient of u 
in the second term on the right-hand side of eq.(10) is positive. 

Katz, et al, deal with this instability by setting this term to zero at 
any point at which this occurs. We deal with it by choosing the 
potential on the outer boundary so that the maximum electron 
relative density obtains on the outer boundary of the calculation 
space. This guarantees that the coefficient, which is the derivative 
of the relative electron density, is negative everywhere. 

Numerical solutions to eqs.(5) have been found in both the 
edge and strip geometries starting with NASCAP/LEO-type 
electrostatic potential solutions resulting from eq.(9) and producing 
the fields, n, F and O' through one iteration cycle. Solutions of 
eq.(9) have been compared with output of the NASCAP/LEO code 
for the electrostatic potential. As expected, the electrostatic 
potentials are quite comparable in regions that are not sensitively 
multi-dimensional. Near the corner in the edge geometry and near 
the voltage discontinuity in the strip geometry, however, 
differences are readily seen. 



3) Boundary Conditions 

Symmetry permits the solutions on the boundaries 
perpendicular to the biased surfaces to be "mirrored”, e.g., 
equipotential lines cross these boundaries perpendicularly. Such 
boundaries include the two side boundaries in the strip geometry 
and the x and y axes in the edge geometry. 

For solutions of Poisson's Equation, the normalized 
electrostatic potential is equal to one on the comer surface in the 
edge geometry (Fig. l.b). The normalized outer boundary potential 
is 0.01, initially. This means that, on a 100 Volt comer, the outer 
boundary would be at 1 Volt. Similarly, the average potential of 
the strip surface (Fig. l.a) is also set to one and the upper boundary 
to 0.01, initially. 

With the outer boundary potential value set initially to 0.01, a 
minimum value of y is found that produces convergence of eq.(9). 
The potential on the outer boundary is then set equal to Ob' - 
(2y) 2 /3. This guarantees that (y - 0' 3 /2/2) £0 everywhere. Solutions 
are then recalculated using this value of Ob’. 

The above process results in an outer boundary potential of 
2-4 Volts for a 100 Volt surface. Thus, since both the ion and 
electron thermal kinetic energies in the ambient plasma are much 
less than 2 e.v., the sheaths modeled in this way are fully Child- 
Langmuir in character and the ions are effectively excluded from 
the sheath. There are no regions in the sheath where eO«kT since 
the electrostatic potential rises monotonically as one proceeds 
inward from the outer sheath to the biased surfaces. 

The boundary conditions on the velocity potential, F, are 
similar to those on the electrostatic potential. The same translation 
and reflection symmetries obtain for F as for O' and so the 
boundaries perpendicular to the biased surfaces are, again, 
"mirrored". However, although definite values of O’ are imposed on 
the biased surfaces and the outer boundary, the values of F on 
these surfaces are not known a-priori. What is known is the 
magnitude of the gradient of F' derived from the Energy Equation, 


(V’F’) 2 = O'. To deal with these boundary conditions, the value of F' 
on the outer boundary is set to zero in our calculations and the 
equation for F' iterated until the Energy Equation is satisfied 
everywhere on the biased surfaces. This requires a delicate 
numerical approach and much computer time. 

4) The Self-Consistent Cycle 

Once F is found using the density function, n, calculated by 
taking the Laplacian of the NASCAP/LEO-type electrostatic 
potential, a new electrostatic potential is calculated from (V’F’) 2 = O'. 
This potential is compared with the NASCAP/LEO-type potential at 
every point in the calculation space. Further processing of the self- 
consistent cycle should consist of testing the agreement of these 
two potentials to some precision and, if they are sufficiently 
different, mixing the potentials, taking the Laplacian of the result so 
as to produce a new relative electron density, n, solving the 
Equation of Continuity for a new velocity potential, F, and finding 
yet another electrostatic potential via the Energy Equation. 

Even after experience with only one processing cycle, it is 
clear that the self-consistent process is fraught with numerical 
instabilities. Solving the Equation of Continuity to high precision is 
especially important. For calculations on a 41x41 point grid, a 
single flow solution can consume many minutes on our VAX 8200 
computer. For this reason, self-consistent solution cycles beyond 
the first cycle have not yet been seriously attempted. However, 
over one cycle it has been found that the input and output 
electrostatic potentials differ significantly near the voltage 
discontinuity in the strip geometry for high-modulation cases and 
almost everywhere in the edge geometry. This means that the 
NASCAP/LEO-type electrostatic potentials are not good solutions of 
the fluid model equations, as suspected. 

5) Comparisons With NASCAP/LEO 

Figures 7&8 show comparisons in the edge geometry between 
an electrostatic potential calculated using NASCAP/LEO and an 
electrostatic potential calculated using the fluid model. We have 
replotted the NASCAP/LEO code output using S. In these figures, 



the surface voltages have been normalized to one. Plasma 
parameters used to calculate these results are as follows: 

Electron Temperature = 0.1 e.v. 

Electron Density = 10 6 /cc 

Debye Length = 0.00235 meters 

Surface Voltage = 100 Volts 

Child-Langmuir Length = 0.525 meters 

As in the Debye Approximation plots (Figs. 2-5), the 
NASCAP/LEO plot of Fig. 7 has a zero on the surface of the 1st outer 
grid. The zero in Fig. 8 is beyond the calculation space, the 
boundary value of the potential having been chosen to maximize 
the electron density on the boundary. One sees that the 
equipotentials in Figs. 7&8 are nearly coincident although the 
NASCAP equipotentials are somewhat compressed toward the 
surface in the outer sheath, no doubt due to the position of the zero 
surface. 

Comparisons of the fluid model and NASCAP/LEO are not 
reported for the strip geometry because of our inability to obtain 
numerically stable runs of NASCAP/LEO using the object shown in 
Fig. 6 with V 00 = 100 Volts. 

6) Fluid Model Runs 

a) Strip Geometry 

Calculations for the strip geometry have been made using the 
plasma parameters given above. Adjacent strips have voltages, V OQ 
plus or minus a percentage modulation. With V Q o — 100 Volts a 
10% modulation means that the higher voltage strips are at 110 
Volts and the lower voltage strips are at 90 Volts. All calculations 
have been done on a 41x41 point grid and, as in the Debye case, 
voltages have been normalized through division by V 0 o* 

Calculations for modulations of 0%, 10%, 20%, 30%, 40%, 50%, 
60%, 70% and 90% have been made. Figures 9-38 show the results 
of calculations having 0%, 10%, 30%, 50%, 70% and 90% modulations. 
The figures are arranged in groups of five. The first figure in each 



group shows the solution of Laplace's Equation. This solution is 
used as the initial electrostatic potential in the iteration of the 
NASCAP/LEO-type, non-linear Poisson's Equation (eq. 9), the 
solution of which is shown in the next figure. 

The third figure shows the normalized electron charge 
density, n = N/N 0 , derived by taking the Laplacian of the Poisson 
solution and multiplying it by the appropriate factor (see eqs. 5). 

The fourth figure is the solution of the Equation of Continuity 
for the velocity potential, F. Fluxes over the outer boundaries and 
onto the biased surfaces have been calculated for these solutions 
and used to check particle conservation. These results are 
summarized in Table I. 

The fifth and last figure in each group shows the output 
electrostatic potential derived from the Energy Equation (also 
denoted as derived from the Velocity Potential). This solution is to 
be compared to the NASCAP/LEO-type, Poisson potential. If these 
potentials are essentially the same, one can conclude that the 
NASCAP/LEO-type potential is a solution of the two-dimensional 
sheath equations presented here. Significant differences between 
these solutions would suggest that the NASCAP/LEO-type potential 
is not a solution of the equations presented. If differences occur 
most prominently near the surface voltage or shape discontinuities, 
one may infer that the one-dimensional charge-potential relation 
assumed in NASCAP/LEO is asserting its presence. 

Figures 9-13 show strip-geometry solutions having no 
modulation. Thus, this case is one-dimensional. Fig. 9 shows a 
Laplace solution that falls-off linearly with distance above the 
uniformly charged surface. This is to be expected since the solution 
of Laplace's Equation in one dimension is a straight line. Fig. 10 
shows the NASCAP/LEO-type potential. One can see that space 
charge pulls the equipotential surfaces toward the biased surface 
somewhat. 

Fig. 11 shows the relative electron density. Because the 
relative electron density is proportional to the Laplacian of the 
electrostatic potential, here it is proportional to the second 



derivative of the potential with respect to y. One can see that the 
second derivative is near zero just above the surface and increases 
as one goes toward increasing y values. 

Fig. 12 shows the velocity potential, F, derived from the 
Equation of Continuity. The velocity field is the gradient of this 
potential. One can see that the velocities are larger near the surface 
than in the outer sheath. This agrees with the behavior of the 
analytical solution of the Equation of Continuity in one dimension, 
namely, nV = const. Where n is small V will be large and where n 
is large V will be small. 

Fig. 13 shows the output electrostatic potential derived from 
the velocity potential. Superposition of Figures 10 and 13 shows 
almost exact agreement between the NASCAP/LEO-type 
electrostatic potential and the new electrostatic potential. Were we 
interested merely in the one-dimensional sheath above a biased 
surface, the self-consistent field process could be considered to 
have converged after only one cycle. This result is to be expected 
since the NASCAP/LEO-type solution should be correct in one 
dimension. 

Figures 14-18 show solutions for a 100 Volt surface with a 
10% modulation. In these figures one can begin to see the effects of 
the differing biases of the strips. Notice that the NASCAP/LEO-type 
Poisson solution is not very different from the Laplace solution. 

The electron densities in the sheath appear so small that, for rough 
purposes, the electrostatic potential may be approximated by the 
Laplace solution. Notice also that the velocity potential shows a 
particle flow only slightly skewed toward the higher voltage strip. 
Apparently, 10% modulation does not change the downward flow of 
electrons appreciably. Finally, notice that the NASCAP/LEO-type 
potential (Fig. 15) and the new electrostatic potential (Fig. 18) 
remain very close. 

Figures 19-38 show results in the strip geometry for a 100 
Volt surface having modulations of 30%, 50%, 70% or 90%. These 
figures show, as expected, that increasing modulation causes 
increased skewing of the potentials toward the higher voltage strip. 
This is also true of the particle fluxes, although the effect is much 



less than would be expected by looking at the velocity potential 
plots and the electron density plots separately. Because flux is the 
product of velocity and density, the high velocities shown over the 
higher voltage strip are almost completely compensated by the low 
densities there. Table I summarizes the associated particle flux 
calculations. Notice in Table I that overall particle loss is extremely 
small in all cases. This means that the Equation of Continuity is 
being solved quite well in these runs. 

Comparisons of the Laplace and the Poisson solutions at each 
value of modulation show that the Laplace solutions are 
increasingly poor approximations to the electrostatic potential in 
the sheath as the modulation is increased. Further, comparison of 
the NASCAP/LEO-type potentials with the potentials derived from 
the velocity potentials shows increasing disagreement as the 
modulation is increased. The assumption presented here is that 
this is due to the breakdown of the one-dimensional charge 
density-electrostatic potential relationship (eq.7) embedded in the 
NASCAP/LEO-type solutions. 

It should be pointed out that, because they differ significantly 
after one cycle of the self-consistent process, neither the 
NASCAP/LEO-type potentials nor the potentials calculated 
subsequently from the velocity potential can be correct solutions of 
the fluid model equations. Correct solutions require convergence of 
the self-consistent field process over, perhaps, many cycles. 
Achieving this awaits a better computing environment. 

b) Edge Geometry 

Figures 39-63 show results of runs made using the edge 
geometry. In all cases the electrostatic potentials are normalized to 
the potential on the surface of the edge, 100 Volts. The voltage on 
the outer boundary is set initially at 1 Volt but is adjusted upward 
to as much as 4 volts for calculational stability. The plasma 
parameters are the same as given in part (B.5) above. 

The figures are arranged in groups of five as above: Laplace 
solution, NASCAP-type solution. Relative Electron Density, Velocity 
Potential and Electrostatic Potential from the Velocity Potential. 



Because there is no easy way to gradually approach the geometry 
of a corner, a series of solutions starting with a one-dimensional 
solution, as was presented for the strip geometry, will not be 
presented. Instead, we present a selection of calculations for 
different sized edges and out to different distances. 

The biased structure shown is 1/4 of a cross section through 
an infinitely long conducting bar. For such an object, a square 
calculation space is inappropriate because the equipotential lines of 
the electrostatic potentials would be pulled artificially toward the 
upper right-hand corner of the calculation space, although they 
would quickly become nearly circular as one proceeded inward. 
Because at large distances one expects the sheath to approximate a 
circle in cross section, anyway, the calculation space has been 
truncated to a circle of 40 grid spaces in radius and centered at the 
point (0,0). Thus the high circularity of potentials of the maximum 
possible radius is also artificial but less so than equipotentials 
calculated using square outer boundaries. 

Figures 39-43 show calculations for 1/4 of a square 
conducting bar 0.346 meters on a side. With Dcl = .5247 meters, 
the calculation space goes out to 1 Dcl from the bar above its flat 
surfaces. The grid is 41x41 points. The voltage on the outer 
boundary is 2.47 Volts. The effective scale length determined from 
the adjusted value of y in the outer sheath is 2 grid spaces. The 
Debye Length is 0.135 grid space so that features of the order of a 
Debye length are unresolved. Table II indicates that good particle 
conservation was achieved for the flow shown in Fig. 42. 

Notice that the NASCAP-type solution (Fig.40) is nearly the 
same as the Laplace solution (Fig.39). However, the electrostatic 
potential derived from the velocity potential (Fig. 43) is pulled very 
much inward relative to the NASCAP-type potential. These two 
potentials are very different, a condition implying that neither is 
the correct solution of the fluid model equations. The correct 
solution must lie in-between these solutions. 

Figures 44-48 show results for a 0.92m x 0.66m bar. The 
calculation space is the same as that of the previous figures. Again 
the Laplace and NASCAP-type solutions are similar although the 



solutions on the lower right boundary are artificially pushed 
toward the bar. In this direction the sheath thickness is forced to 
be less than 1 Dcl- 

As in the previous case, particle conservation is obtained to 
high precision. Also, the ratio of fluxes to the top and to the side of 
the corner are very nearly the same as the ratio of the extents of 
these surfaces, i.e., 14/10. As before, the electrostatic potential 
derived from the velocity potential is much pulled-in relative to the 
NASCAP-type potential. 

Figures 49-53 show the same results as do Figures 39-43 but 
a different scale is used. The grid spacing is 0.05 Dcl in these 
figures rather than 0.033 Dcl- Thus the bar's dimensions are 0.525 
m on a side. The calculation space goes out to 1.5 Dcl = 0.787 m 
above the flat surfaces of the bar. For calculational stability, the 
outer boundary has been set to 4.3 Volts. 

Notice in the first two of these figures that the Laplace and 
NASCAP-type solutions are no longer similar. As expected, the 
Laplace solution extends out to the calculational boundary but the 
NASCAP-type solution of Fig. 50 terminates before this, at about 
0.95 Dcl- The sheath radius compares to a radius of 0.86 Dcl in the 
NASCAP-type solution of Fig. 40. The slightly greater sheath radius 
in Fig. 50 is perhaps due to the larger boundary voltage, which 
would force the equipotentials outward somewhat. Nevertheless, 
the calculations shown in Fig. 50 have resulted in a sheath of the 
expected, finite thickness, namely, about 1 Dcl- 

Reference to Table II shows that the flow of Fig. 52 conserves 
particles very well. Fig. 53 shows, again, that the final electrostatic 
potential and the NASCAP-type potential are very different. The 
wiggles on the outer equipotential in Figures . 53 and 49 do not 
represent real variation. They are probably due either to plotting 
artifacts, the fact that a smooth, circular boundary cannot be 
achieved on a square grid of points or budding numerical 
instabilities in the outer sheath. 

Figures 54-58 show the same kind of information as do 
Figures 44-48 but with a grid spacing of 0.05 Dcl. Thus the bar size 



in these figures is 0.735 m x 0.525 m. The calculation space is the 
same as that of Figures 44-48. The outer boundary voltage is 4.29 
Volts. 


One sees in these figures that the NASCAP-type solution 
shows a sheath of definite radius, about 1 Dcl. and that the final 
potential disagrees strongly with the NASCAP-type potential. As 
before, the Equation of Continuity is well solved and particles are 
conserved to high precision. Note, particularly, that in Fig. 55 the 
outer equipotentials are very circular even though the calculation 
space boundary is quite a bit further out. Indeed, this circularity is 
not lost until one approaches to within about 0.5 Dcl of the corner. 

Figures 59-63 show calculations for a square bar of the same 
dimensions as that of the bar shown in Figures 39-43, i.e., 0.66 Dcl 
on a side. In these figures, however, the calculation space is 
extended to 1.56 Dcl above the flat surfaces of the bar. For 
numerical stability, the effective scale length of the outer sheath 
has been increased from 2 to 2.1 grid spaces. This means that the 
outer boundary electrostatic potential is 4.2 Volts. 

One sees in the NASCAP-type solution of Fig. 60 that the 
sheath extends out only to about 0.71 Dcl above the flat surfaces. 
This is to be compared to Fig. 40, in which the sheath extends out to 
about 0.86 Dcl but the calculation space extends only out to 1 Dcl- 
This difference may be due to proximity of the calculation 
boundary to the outer sheath surface in Fig. 40, i.e., the boundary 
may be pulling the sheath outward. The sheath shown in Fig. 60 is 
less affected by the boundary, which is farther away. One should 
note here once again that the boundary of the calculation space is a 
circle, 40 grid spaces in radius and centered at (0,0). 

As in the other cases presented, Fig. 62 shows a flow that 
conserves particles quite well (see Table II). Also, the final 
potential is very much different from the NASCAP-type solution, as 
before. 


C. Summary and Conclusions 

The modeling of the Debye Approximation electron sheaths 
in the edge and strip geometries was completed .in the summer of 
1989. Electrostatic potentials in these sheaths were compared to 
NASCAP/LEO solutions for similar geometries. Velocity fields, 
charge densities and particle fluxes to the biased surfaces were 
calculated for all cases. 

The major conclusion to be drawn from the comparisons of 
our Debye Approximation calculations with NASCAP-LEO output is 
that, where comparable biased structures can be defined and 
sufficient resolution obtained, these results are in general 
agreement. 

Numerical models for the Child-Langmuir, high-voltage 
electron sheaths in the edge and strip geometries were constructed 
tn-4989. Electrostatic potentials were calculated for several cases in 
each of both geometries. Velocity fields and particle fluxes were 
calculated. The self-consistent solution process was carried through 
one cycle and output electrostatic potentials compared to NASCAP- 
type input potentials. 

The major conclusions to be drawn from the Child-Langmuir 
analysis of the strip and edge geometries are as follows: 

1) Equations for the electron sheath that connect the electrostatic 
potential, electron density and velocity field can be derived from a 
multi-dimensional, warm fluid plasma model. These equations 
reduce to those used in the NASCAP-LEO code in one dimension. 

2) NASCAP-type electrostatic potentials, which are solutions of the 
multi-dimensional Poisson’s Equation but using the one- 
dimensional relation between charge density, velocity and 
electrostatic potential, are not compatible with the warm plasma 
fluid model near points of severe voltage or shape discontinuity. 

The incompatibility becomes more pronounced as voltage or shape 
discontinuities are made larger. 



3) If particle velocity fields are assumed irrotational, the Equation 
of Continuity becomes an equation for the velocity potential, which 
can be solved to high precision. This requires imposing a non- 
linear boundary condition on the biased surface. The Continuity 
Equation approach to determining particle flow is an alternative to 
particle tracking and has the feature that it embodies particle 
interactions during the flow. 

4) Sheath thicknesses of the NASCAP-type solutions in the edge 
geometry are approximately equal to but somewhat less than the 
expected value of one Child-Langmuir length. However, the 
electrostatic potentials calculated from the velocity fields have 
much smaller sheath thicknesses. If the self-consistent solution 
cycle were to be pursued to convergence, one might expect that the 
final solutions would show intermediate sheath thicknesses. This 
would imply generally thinner sheaths than are predicted by 
NASCAP/LEO-type solutions. 
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Table I 

Particle Fluxes for the Strip Geometry 

100 Volt Surface +/- Modulation 
(arbitrary units) 


Fig.# 

% Mod. 

LFlux 

RFlux 

TFlux 

%NC 

%RtShift 

12 

0 

.1783 

.1783 

.3565 

-.03 

0 

17 

1 0 

.1763 

.1757 

.3519 

-.03 

-.2 

22 

30 

.1748 

.1756 

.3503 

-.03 

.2 

27 

50 

.1722 

.1745 

.3467 

0 

.6 

32 

70 

.1669 

.1722 

.3391 

0 

1.6 

37 

90 

.1514 

.1671 

.3187 

+.06 

5.0 


LFlux = integrated normal flux on lower left (low voltage) boundary 

RFlux = integrated normal flux on lower right (high voltage) boundary 

TFlux = integrated normal flux on upper boundary 

%NC = 100*[l-(lflux+rflux)/topflux] = % non-conservation of flux 

%RtShift = 100*(rflux-lflux)/topflux = % of flux shifted to higher voltage 
surface 
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Table H 

Particle Fluxes for the Edge Geometry 
(arbitrary units) 


Fig.# 

CTFlux 

CCFlux 

CRFlux 

TTFlux 

TRFlux 

%NC 

42 

.084 

.0063 

.084 

.088 

.088 

.5 

47 

.117 

.0063 

.085 

.104 

.105 

.3 

52 

.084 

.0063 

.084 

.088 

.088 

.5 

57 

.117 

.0063 

.085 

.104 

.105 

.4 

62 

.059 

.0063 

.059 

.063 

.063 

.8 


CTFlux = integrated normal flux onto top of edge 
CCFlux = flux to comer 

CRFlux = integrated normal flux onto right side of edge 

TTFlux = integrated normal flux crossing top boundary 

TRFlux = integrated normal flux crossing right boundary 

%NC = 100*[1 -(total outer boundary flux/total flux to edge) = % non- 
conservation of flux 
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